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1.0  INTRODUCTION 

Turbulent  ducted  flow  is  commonly  found  in  diffusers,  nozzles,  test  cells,  and  wind 
tunnels.  The  complexity  of  the  flow  involved  depends  upon  many  factors  such  as  the 
geometry  of  the  wall,  the  turbulent  nature  of  the  flow,  flow  separation,  three-dimensional 
effects,  the  Mach  number  effect,  etc. 

For  incompressible  turbulent  flow,  a method  (Refs.  1 and  2)  has  been  developed  to 
obtain  the  finite-difference  solution  of  the  Navier-Stokes  equations  with  a two-equation  k-e 
turbulence  model.  It  has  been  used  in  the  computation  of  both  separated  and  nonseparated 
turbulent  diffuser  flows  with  nonuniform  inlet  conditions.  The  original  analysis  was 
formulated  for  two-dimensional  (2-D)  or  axisymmetrical  flows  with  a single-wall  boundary 
layer  and  a coordinate  transformation.  In  order  to  handle  general  nonsymmetrical  2-D  flows 
and  annular  diffuser  flows,  another  wall  boundary  layer  and  coordinate  transformation 
were  added  to  the  analysis  in  the  present  study.  With  this  modification,  a fairly  general  class 
of  ducted  flow  problems  can  now  be  computed  easily  with  the  incompressible  Navier-Stokes 
code.  Among  other  important  changes  are  (1)  a simple  exponential  coordinate 
transformation  to  handle  the  two  boundary  layers,  (2)  a formula  for  the  initial  flow-field 
guess,  and  (3)  a new  approach  for  computing  the  pressure  field. 

For  the  compressible  flow  computations,  the  requirement  to  calculate  the  density 
variation  makes  it  necessary  to  know  the  temperature  and  pressure  fields.  It  is  the  purpose  of 
the  present  study  to  extend  the  incompressible  analysis  to  the  compressible  flow  regime 
within  the  framework  of  the  vorticity-stream  function  formulation  of  the  Navier-Stokes 
equations.  The  derivation  of  the  governing  equations  in  terms  of  the  vorticity  and  the  stream 
function,  the  coordinate  transformations,  the  logic  of  the  numerical  procedure,  and  some 
preliminary  results  are  presented. 

2.0  THE  VORTICITY-STREAM  FUNCTION  FORMULATION 

The  basic  governing  equations  are  the  continuity,  momentum,  and  energy  equations. 
The  set  of  steady-state  equations  can  also  be  expressed  in  terms  of  the  vorticity  and  the 
stream  function.  These  are  given  in  both  2-D  Cartesian  and  cylindrical  coordinates  as 
follows  (different  vorticity-stream  function  formulations  can  be  found  in  Refs.  3 and  4): 

Vorticity  Equation 
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Stream  Function  Equation 
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In  the  vorticity  equation,  the  first  two  second-order  terms  represent  the  diffusion  term,  the 
next  two  first-order  terms  represent  the  convection  term,  and  the  rest  represent  the  vorticity 
source  terms.  Of  the  vorticity  source  terms,  the  first  term  comes  from  the  cylindrical  coor- 
dinates used,  the  second  term  is  related  to  the  density  variation,  the  next  two  terms  are 
related  to  the  first-order  eddy  viscosity  variation,  and  the  last  term  is  related  to  the  second- 
order  eddy  viscosity  variation.  The  eddy  viscosity  (^)  and  the  density  (p)  which  appear  in  the 
vorticity  and  the  stream  function  equations  must  be  determined  from  the  turbulence  model- 
ing, the  energy  equation,  the  equation  of  state,  and  the  pressure  equation. 

3.0  THE  PRESSURE  EQUATION 


There  are  several  methods  available  to  recover  the  pressure  field,  namely  (1)  by  direct 
integration  of  the  momentum  equations,  (2)  by  solving  a second-order  pressure  equation, 
and  (3)  by  solving  a new  pressure  equation  described  in  the  present  analysis. 
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3.1  The  Conventional  Pressure  Equation 

i In  the  first  approach,  it  is  possible  to  obtain  a pressure  field  distribution  by  direct 
integration  of  the  momentum  equation.  In  the  region  where  the  flow  variables  are  relatively 
smooth  and  continuous,  the  pressure  integration  is  straightforward.  Otherwise,  the  local 
error  generated  will  be  carried  along  the  path  of  integration.  As  a result,  the  pressure  level  at 
the  end  of  the  integration  path  is  usually  inaccurate.  This  happens  often  along  the  wall 
where  the  velocity  gradient  is  large.  This  method  can  be  used  when  the  velocity  field 
obtained  is  accurate  and  there  are  enough  grid  points  to  define  the  solution  smoothly.  The 
line  integration  methods  used  to  obtain  the  pressure  distribution,  although  good  for  some 
simple  applications,  do  not  generally  provide  a unique  solution. 

In  the  second  approach,  a Poisson  type  of  pressure  equation  is  derived  from  the 
momentum  equation  through  a (V  *)  operation.  The  result  is  a second-order  field  equation, 

( 52p  <32p\  § / <?p\ 

\d*2+  <9r2/  r \dr)  ~ ‘P  (3> 

where 


For  2-D  incompressible,  constant  viscosity  flow,  Eq.  (1)  is  reduced  to  the  following  form: 
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where  p is  defined  as  (p/p). 

The  second-order  pressure  equation  derived  represents  an  improvement  over  the  direct 
integration  scheme.  But  there  are  difficulties  associated  with  the  boundary  conditions.  In 
most  ducted  flow  problems,  the  pressure  along  the  wall  is  not  known.  For  flow  between  two 
walls,  the  zero  normal  pressure  gradient  condition  would  have  to  be  used.  As  a result,  the 
boundary  condition  becomes  weak  in  controlling  the  absolute  pressure  level.  In  addition, 
slow  convergence  characteristics  usually  exist  for  the  gradient-type  boundary  conditions.  So 
far,  no  satisfactory  method  has  been  developed  to  improve  the  situation. 

3.2  A New  Pressure  Equation 

With  the  weakness  of  the  second-order  pressure  equation  known,  it  is  possible  to  create  a 
new  pressure  equation  from  the  momentum  equations  with  stronger  convergence  properties. 
The  idea  is  to  retain  as  much  as  possible  the  first-order  nature  of  the  momentum  equations, 
while  at  the  same  time  working  with  a second-order  field  equation.  Although  the  approach  is 
highly  heuristic,  the  results  of  this  approach  discussed  in  later  sections  do  support  the  basic 
idea.  The  basic  idea  behind  this  new  equation  is  simple  and  is  explained  briefly  as  follows.  In 
most  ducted  flow  problems,  there  are  inviscid  and  viscous  boundary-layer  regions.  In  the 
subsonic  inviscid  region,  the  flow  variables  change  relatively  smoothly  along  the  main  flow 
direction  {or  along  the  streamline).  It  is  possible  to  determine  a path  function  in  such  a way 
that  the  resultant  pressure  equation  is  highly  convective  along  this  path.  With  proper 
selection  of  the  coefficient,  one  can  make  the  second-order  diffusion  term  relatively  small. 
Thus  one  can  achieve  the  forward  integration  of  the  pressure  by  solving  a second-order 
diffusion-convection  pressure  equation.  In  order  to  obtain  the  pressure  field  for  the  rest  of 
the  flow  field,  one  can  specify  the  path  function  normal  to  the  main  flow  direction  (or 
normal  to  the  boundary  layer). 

The  pressure  equation  with  this  new  set  of  path  functions  will  pick  up  the  pressure 
information  along  the  main  flow  direction  and  propagate  it  across  the  boundary  layer  to 
reach  the  wall.  This  is  a natural  approach,  in  a sense,  because  the  normal  pressure  gradient  is 
usually  small.  Therefore,  one  would  expect  to  have  relatively  accurate  pressure  distributions. 
The  key  to  the  success  is  the  creation  and  specification  of  the  vector  path  function.  The  new 
pressure  equation  derived  from  the  momentum  equations  with  the  vector  path  functions  PI 
and  P2  is  as  follows: 

\ dr2  J r \<3r  / {fi  J-  \dx  / ^ Si  J 
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where  PI  and  P2  are  the  vector  path  functions  and  Ug  is  a reference  velocity. 


The  strong  convective  nature  of  this  new  second-order  pressure  equation  is  controlled  by 
the  path  function  and  the  coefficient  [puc/(/i  + The  path  function  must  be  specified  in 
advance.  It  could  be  changed  to  improve  the  accuracy  on  the  basis  of  the  velocity  dis- 
tribution. In  general,  the  path  function  depends  on  a particular  flow  problem  to  be  solved. 
For  the  ducted  flow  application,  PI  is  one  along  the  centerline  or  the  midposition  line  and  is 
zero  elsewhere.  This  has  an  effect  of  convecting  the  pressure  along  the  midposition  line  in 
the  downstream  direction.  The  other  path  function,  P2,  is  zero  along  the  midposition  line,  is 
one  in  the  positive  half  of  the  flow  field,  and  is  minus  one  in  the  lower  half  of  the  region. 
The  information  about  the  pressure  is  then  spread  away  from  the  midposition  line  in  the 
normal  direction.  The  underlying  physics  is  explained  in  the  following  example: 
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For  a parallel  channel  flow  with  a boundary  layer  and  a core  region,  the  new  pressure 
equation  behaves  in  the  core  region  as 


(8 


(6) 


On  the  other  hand,  in  the  boundary  layer  region,  it  becomes 


f^p\ 

\<W 


. s„ 

dr  1 P 

^ / 


(7) 


Because  of  the  strong  coefficient  (puc/»,  the  pressure  field  is  being  convected  along  the 
prescribed  direction  of  PI  and  P2. 

Several  obvious  advantages  can  be  obtained  by  using  the  newly  derived  pressure 
equation:  (1)  the  second-order  diffusion-convection  equation  can  be  easily  solved  by  the 
general  finite-difference  formulation  with  decay  function;  (2)  the  path  function  can  be 
selected  in  such  a way  that  the  error  in  the  pressure  field  is  minimized;  (3)  a fast  convergence 
rate  is  realized  because  of  the  strong  convective  nature  of  the  equation;  and  (4)  the  problem 
with  the  wall  boundary  condition  is  relieved  because  the  wall  pressure  is  now  approached  in 
the  normal-to-the-wall  direction. 


Some  feeling  for  the  complexity  of  the  transformed  equation  in  the  computational 
domain  may  be  gained  from  the  transformed  pressure  equation  shown  here. 


Transformed  Pressure  Equation  for  Compressible  Flow 
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4,0  A TWO-EQUATION  k-e  TURBULENCE  MODEL 


(8) 


The  closure  problem  related  to  the  turbulent  shear  stresses  is  handled  in  the  present 
analysis  through  the  eddy  viscosity  concept.  The  eddy  viscosity  (fit)  is  determined  from  a 
two-equation  k-e  model  through  the  Prandtl-Kolmogorov  relation.  The  molecular  viscosity 
effect,  which  is  important  in  the  viscous  sublayer  region,  is  modeled  through  the  eddy 
viscosity  coefficient  (C^)  and  additional  terms  in  turbulent  kinetic  energy  (TKE)  and 
dissipation  equations. 

The  Prandtl-Kolmogorov  Relation  for  #it  is 

k2 

^t-^7 


where 


C = 

**  3(a  +•  A/b) 


, A = \J  2 Icy  pi  ft,  a = 1,100,  b = 0.27 


(9) 
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and  where  k is  determined  from  the  turbulent  kinetic  energy  equation: 

i L -*i).* — 1 — 

\(?x2  <?r2/  (p  + pt)  \ dx  } <9x  (ft  + fit)  j_  dr  r J dr 

(10) 


+ |_\d«/  \dr  / \ r / \dr  dx/  ( (ft  + pt)\ 


and  e is  determined  from  the  TKE  dissipation  equation: 


(11) 

where 

Ct  = 1.36,  C2  = 1.92  [l-0.3  exp  <-R2)],  R = pk2/p<). 

5.0  BOUNDARY  CONDITIONS  AND  INLET  CONDITIONS 

The  boundary  conditions  for  turbulent  ducted  flow  problems  can  be  divided  into  four 
groups,  namely  wall  condition,  symmetry  condition,  exit  condition,  and  inlet  condition.  Of 
these,  the  inlet  condition  requires  special  attention.  Therefore,  it  is  discussed  separately 
below.  For  the  solid-wall  boundary  condition,  the  nonslip  condition  is  usually  used  (i.e.,  u 
= v = 0).  When  a line  of  symmetry  exists  in  the  computational  domain,  a zero  normal 
gradient  condition  is  commonly  used.  At  the  downstream  side  of  the  computational  domain, 
it  is  necessary  to  specify  the  exit  conditions  because  of  the  elliptic  nature  of  the  problem. 
Since  one  does  not  know  the  exit  solution  in  advance,  a zero  axial  gradient  condition  is 
normally  used  to  simulate  the  nearly  parallel  condition  at  the  exit  station. 

Inlet  conditions  for  the  turbulent  ducted  flow  problem  normally  consist  of  a core  region 
and  boundary-layer  regions.  The  profiles  in  these  regions  usually  are  highly  nonuniform. 
They  can  be  either  obtained  from  experiment  or  specified  by  analytical  modeling.  In  most 
cases,  it  is  not  possible  to  obtain  the  complete  set  of  data  needed  for  computational  purposes 
from  the  core  region  into  the  turbulent  boundary  layer-region  and  through  the  relatively  thin 
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sublayer  to  the  wall.  Therefore,  it  is  highly  desirable  to  have  formulae  which  can  be  used  to 
generate  a set  of  complete  profiles  for  the  numerical  computation  from  some  gross  features 
of  the  inlet  flow.  In  the  present  analysis,  the  inlet  condition  is  analytically  calculated  on  the 
basis  of  the  given  boundary-layer  thicknesses  and  the  Reynolds  number.  Within  the 
boundary-layer  region,  the  modified  Van  Driest  formula  is  used  to  provide  the  complete 
profile  of  the  velocity  gradient  throughout  the  viscous  layer  on  the  basis  of  an  assumed  total 
shear  stress  distribution  (Ref.  2).  The  related  flow  variables  such  as  velocity  and  stream 
function  are  derived  numerically. 

With  a two-equation  turbulence  model,  it  is  necessary  to  know  the  turbulent  kinetic 
energy  (k)  and  the  dissipation  rate  (e)  throughout  the  whole  turbulent  boundary  layer.  This 
cannot  be  easily  achieved  experimentally.  Therefore,  it  is  computed  analytically  from  the 
following  TKE-shear  stress  relation 


k = 

p \dy  ) 


(12) 


where  is  modelled  as  A/3  (a  + A/b),  A = a/2E  y/v,  a = 1,100  and  b = 0.27.  The  above 
equation  is  an  equation  for  k and  y.  It  must  be  solved  iteratively  because  of  its  nonlinearity. 
The  asymptotic  behavior  of  this  relation  in  the  region  away  from  the  thin  wall  layer  is  given 
by  the  conventional  TKE-shear  stress  formula  as 

k . _L 

0.3  p \dyf  (13) 

It  must  be  emphasized  again  that  it  is  very  important  to  provide  a self-consistent  set  of 
profiles  so  that  one  can  avoid  sudden  readjustments  of  the  numerical  solution  near  the  inlet 
station. 


6.0  NUMERICAL  METHOD  AND  COORDINATE  TRANSFORMATION 

In  the  present  analysis,  the  physical  flow  field  is  first  transformed  into  a rectangular 
computational  domain,  and  the  solution  is  obtained  iteratively,  with  an  initial  flow-field  ' 
guess,  through  a general  finite-difference  formulation  with  decay  functions. 

6.1  COORDINATE  TRANFORMATION 

The  general  coordinate  transformation  is  described  in  Refs.  1 and  2.  In  short,  the  domain 
of  interest  is  first  transformed  into  a rectangular  region  through  a nonorthogonal  coordinate 
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transformation.  A second  coordinate  stretching  is  then  applied  to  provide  an  adequate 
solution  definition  throughout  the  boundary  layer,  including  the  sublayer.  In  the  former 
approach,  the  coordinate  stretching  was  provided  by  two  matched  stretching  functions, 
while  in  the  present  approach,  a simple  exponential  coordinate  stretching  is  employed  for 
the  boundary-layer  region.  The  difference  is  in  the  ease  of  handling.  In  the  core  region,  a 
simple  uniform  grid  arrangement  is  used. 

K * 1 

y = AN (in  the  boundary-layer  region) 

K - 1 

Ay  = Constant  (in  the  core  region) 

where  y is  measured  normal  to  the  wall,  Ah|  is  the  first  mesh  size,  and  K is  the  ratio  between 
two  successive  mesh  sizes.  The  first  mesh  size,  Ah1(  can  be  determined  by  selecting  the  first 
grid  point  inside  the  viscous  sublayer  at  y+  =2  (i.e.,  hi  = 2 p/v*).  The  coefficient  K must 
be  determined  iteratively  when  the  range  of  y and  the  total  number  of  grid  points  in  the 
boundary-layer  region  are  given. 

For  example,  the  condition  is 


where  JN  is  the  number  of  grid  points 
transformation  factors  are  evaluated  numerically  after  the  transformation  function  is 
determined.  The  transformation  function  has  the  advantage  of  smoothing  out  the 
discontinuity  which  normally  occurs  at  the  matching  points.  It  must  be  emphasized  that  the 
proper  evaluation  of  the  transformation  factors  is  important  and  that  the  reproduction  of 
the  inlet  profile  can  be  used  to  check  the  validity  of  the  transformation. 

6.2  INITIAL  FLOW-FIELD  GUESS 

The  initial  flow- field  guess  for  the  steady-state  iterative  procedure  is  very  important.  It 
affects  not  only  the  rate  of  convergence  but  also  the  stability  of  the  iteration  procedure.  For 
turbulent  flow  computations  with  a low  Reynolds  number  model,  it  is  not  known  whether 
there  is  a unique  path  from  an  initially  laminar  region  to  a finally  turbulent  region,  and  vice 
versa.  Therefore,  it  is  important  to  use  an  initial  guess  with  turbulence  properties  which  are 
reasonably  close  to  the  final  solution. 


2v  KJN  - 1 


v*  K — 1 


(15) 


which  covers  the  range  yma*.  The  coordinate 
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In  the  present  approach,  the  inlet  condition  is  used  to  obtain  the  initial  flow-field  guess. 
The  local  profile  is  decomposed  into  two  parts,  namely  the  inlet  profile  part  and  a wake 
profile  part.  These  two  profiles  are  combined  to  satisfy  the  global  continuity  equation  (Fig. 
1).  For  example,  the  velocity  profile  is  written  as 


At  the  inlet,  Au  = 0 and  umean  = ui,c;  therefore,  u = uKy). 


(16) 


Figure  1.  Decomposition  of  initial  velocity  profiles. 


With  the  velocity  profile  decomposition  method,  the  initial  stream  function  distribution 
is  derived.  The  related  variables  such  as  the  vorticity  and  the  velocity  are  obtained 
numerically.  The  initial  viscosity  distribution  is  computed  as 


-Au)  . 

, „ . \JUSJU! 1 + 0.04  • Ay.  ■ Au  • 

1 1 » 1 11  _ 1 


(17) 


where  the  first  term  represents  the  contribution  from  the  adjusted  inlet  condition  and  the 
second  term  represents  the  corresponding  wake  contribution;  Ay[  is  the  assumed  local 
boundary-layer  thickness.  Improvement  in  the  initial  flow-field  guess  could  be  further 
carried  out  by  combining  experimental  data,  computed  results,  and  a simple  integral  method 
to  make  the  Navier-Stokes  computer  code  “smart.” 
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6.3  NUMERICAL  SOLUTION  PROCEDURE 

The  governing  equations  and  the  corresponding  boundary  conditions  are  solved  in  the 
transformed  space  by  a general  finite-difference  method  with  decay  functions  (Refs.  5 and 
6).  The  program  logic  for  the  incompressible  flow  computation  is  shown  in  Fig.  2.  The 
vorticity  and  the  stream  function  equation  form  a main  loop  along  with  the  eddy  viscosity 
computation.  The  pressure  equation  forms  another  separate  loop  because  of  the  pressure 
decoupling  mechanism  in  incompressible  flow. 

For  the  compressible  flow  computation,  the  vorticity,  stream  function,  eddy  viscosity, 
energy,  and  pressure  equations  form  a main  loop  (Fig.  3).  The  pressure  field  must  be 
updated  constantly  because  of  the  density  coupling  in  the  compressible  flow.  Basically,  the 
compressibility  effect  appears  in  the  convection  term  because  of  the  density  modification. 

7.0  RESULTS  AND  DISCUSSION 

The  application  of  the  general  finite-difference  method  with  decay  functions  to  the 
solution  of  the  Navier-Stokes  equations  in  terms  of  the  vorticity-stream  function 
formulation  is  new.  Therefore,  it  is  very  important  to  know  the  validity  and  the  limitation  of 
this  approach.  In  previous  work  (Refs.  1 and  2)  it  was  studied  very  carefully  and  was  applied 
to  the  laminar  diffuser  flow  computations,  where  the  numerical  solution  can  be  easily 
verified  without  the  complication  of  turbulence  modeling.  With  the  coordinate  stretching, 
the.  technique  was  applied  to  the  computation  of  incompressible  2-D  planar  diffuser  flows 
and  several  separated  and  nonseparated  conical  diffuser  flows.  The  results,  by  comparison 
with  a limited  amount  of  available  data,  were  encouraging.  Obvious  extensions  to  the 
previous  work  were  1)  to  include  more  complicated  boundaries,  such  as  annular  diffusers, 
and  2)  to  account  for  compressibility  effects.  To  evaluate  the  first  extension,  the  analysis  and 
the  computer  code  were  applied  to  a wind  tunnel  diffuser  with  centerbody.  To  answer  the 
question  about  compressibility,  several  steps  were  taken.  These  are  described  in  the 
following  sections. 

7.1  INCOMPRESSIBLE  WIND  TUNNEL  DIFFUSER  FLOW  COMPUTATION 

The  wind  tunnel  diffuser  modeled  is  part  of  the  AEDC  Propulsion  Wind  Tunnel  Facility 
(PWT)  Aerodynamic  Wind  Tunnel  (16T).  The  outer  wall  of  the  diffuser  is  three- 
dimensional.  A square  inlet  section  is  followed  by  a transition  section  which  changes  the 
cross  section  from  square  to  circular  (Fig.  4).  The  diffuser  centerbody  is  axisymmetric.  The 
computations  were  made  for  an  equivalent  axisymmetric  annular  diffuser  which  has  the 
same  axial  distribution  of  flow  area  as  the  actual  diffuser.  The  computed  flow  field  can  be 
expected  to  be  significantly  in  error  in  those  parts  of  the  diffuser  where  the  cross  section  of 
the  outer  wall  deviates  greatly  from  axisymmetric.  However,  if  the  predicted  flow  field  is  in 
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Figure  2.  Flow  chart  for  incompressible 
flow  solution  procedure. 


Data  and  Control  Parameters 


Figure  3.  Flow  chart  for  compressible 
flow  solution  procedure. 
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Figure  4.  16T  diffuser  geometry. 

reasonable  agreement  with  experiment,  then  the  numerical  computations  can  be  used  to 
assess  the  effect  of  geometric  and  Reynolds  number  variations  on  the  performance  of  the 
16T  diffuser. 

In  the  course  of  this  analysis  a completely  new  formulation  was  carried  out  to  include 
both  the  wall  and  the  centerbody  turbulent  boundary  layers.  Unfortunately,  this  placed  a 
very  heavy  demand  on  the  coordinate  transformation,  because  it  meant  that  one  must 
handle  two  viscous  layers,  two  sublayers,  and  a core  region.  The  two-function  matching 
procedure  proved  to  be  insufficient  to  provide  a smooth  and  continuous  transformation  in 
the  realistic  Reynolds  number  range"  of  4 x 106  with  about  61  grid  points  in  the  radial 
direction.  This  problem  was  later  solved  with  a simple  exponential  function  transformation 
which  is  described  in  Section  5.1.  With  this  approach,  the  core  region  was  handled  by  a 
uniform  grid,  and  the  two  boundary  layers  were  covered  by  a two-exponential  function 
transformation. 

In  the  computation,  a 69  by  105  grid  was  generated  in  the  radial  and  axial  directions.  In 
the  radial  direction,  the  two  boundary  layers  along  the  outer  diffuser  wall  and  the 
centerbody  were  covered  by  the  grid.  The  centerbody  was  extended  to  the  inlet  plane  with  a 
small  sting  to  eliminate  the  direct  computation  of  the  stagnation  point  region  of  the 
centerbody.  The  computed  wall  and  the  centerbody  pressure  coefficients  are  shown  in  Fig. 
5.  Also  in  Fig.  5,  experimental  results  are  presented  from  a 1/16  scale  model  of  the  diffuser 
tested  in  the  AEDC-PWT  Aerodynamic  Wind  Tunnel  (IT).  Agreement  with  IT  model  data 
is  good  except  for  the  wall  pressure  in  the  square  inlet  section,  where  large  three-dimensional 
effects  exist.  The  Reynolds  number  per  foot  is  4 x 106.  The  computations  required 
approximately  one  hour  on  an  IBM  370/165  computer.  The  storage  requirement  with 
double  precision  is  about  1.2  x 106  bytes.  From  the  convergence  history  shown  in  Fig.  6,  it 
can  be  seen  that  the  solution  reached  a converged  state  in  about  30  minutes.  The  pressure 
distribution  was  obtained  from  the  new  pressure  equation.  The  convergence  history  with  the 
new  equation  is  shown  in  Fig.  7.  It  can  be  seen  that  the  convergence  is  relatively  fast 
compared  to  the  vorticity-stream  function  cycle.  This  is  very  important  because  the  pressure 
cycle  must  be  updated  constantly  in  the  future  computation  of  compressible  flows. 
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b.  Scoop  (centerbody)  Cp  distribution 
Figure  5.  Cp  distribution  for  16T  model  diffuser. 
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Figure  7.  Convergence  characteristics  for  the  new  pressure  equation. 
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Obviously  the  ability  to  efficiently  obtain  accurate  pressure  distributions  with  the  new 
method  has  been  successfully  demonstrated.  Finally,  it  is  worthwhile  to  mention  that  these 
computations  were  carried  out  with  a uniform  axial  mesh  arrangement,  which  may  not  be 
adequate  for  flows  with  large  variation  in  geometry  along  the  axial  direction.  Therefore, 
extension  to  include  a variable  mesh  in  the  axial  direction  is  also  desirable. 

Although  only  the  pressure  distributions  are  presented  here,  a complete  set  of  flow-field 
variables  was  tabulated  as  the  output  of  the  computer  code. 


7.2  COMPRESSIBLE  LAMINAR  FLOW  IN  A SUDDEN  EXPANSION  DIFFUSER 

Although  the  vorticity-stream  function  formulation  has  been  very  successful  in  handling 
incompressible  two-dimensional  flow  problems,  the  capability  of  this  approach  to  obtain 
compressible  flow  solutions  must  be  demonstrated.  The  reason  behind  this  is  that  it  becomes 
necessary  to  update  the  pressure  constantly  by  solving  a pressure  equation.  The  old  problem 
associated  with  the  conventional  pressure  equation  (i.e.,  slow  convergence  and  boundary 
condition)  has  discouraged  people  from  using  it.  It  is  therefore  very  important  to  select  a 
well-defined  and  physically  realistic  test  case  so  that  the  present  logic  for  the  compressible 
flow  can  be  verified.  For  this  reason,  the  planar  laminar  flow  in  a sudden  expansion  diffuser 
was  chosen  to  eliminate  the  additional  complication  involved  in  the  use  of  the  k-e  turbulence 
model. 

The  velocity  profile  is  shown  in  Fig.  8.  At  the  inlet,  the  Mach  number  in  the  core  region  is 
0.8,  and  there  is  a laminar  boundary  layer  near  the  wall.  The  computation  was  started  at  one 
step  height  (H)  ahead  of  the  expansion  so  that  the  elliptic  nature  of  the  flow  could  be 
properly  included.  The  flow  separates  smoothly  from  the  corner  and  forms  a long  separation 
region.  For  a Reynolds  number  equal  to  70,  the  length  of  the  recirculation  region  is  about 
7.5H.  The  centerline  velocity  also  decreases  continuously  because  of  the  diffusion  process. 

In  the  computation,  the  conventional  pressure  equation  of  the  Poisson  type  was  used.  It 
is  likely  that  one  can  include  the  new  pressure  equation  in  the  analysis  and  modify  the 
computer  code  to  obtain  accurate  pressure  distributions  in  a more  efficient  way.  The 
convergence  characteristics  are  shown  in  Fig.  9 with  the  conventional  pressure  equation.  The 
grid  arrangement  is  a 101  by  36  uniform  mesh,  and  the  computing  time  was  about  8.5  min. 
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U 


Axial  Distance,  x/H 


Figure  8.  Velocity  and  streamline  distribution  of  a sudden  expansion  diffuser. 


Figure  9.  Convergence  characteristics  of  sudden  expansion  diffuser  computation. 
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8.0  CONCLUDING  REMARKS 

The  present  steady-state  vorticity-stream  function  formulation  has  been  shown  to 
provide  convergent  solutions  for  both  compressible  and  incompressible  internal  flows.  The 
new  pressure  equation  derived  also  provides  fast  and  accurate  pressure  distributions  for  the 
fairly  complicated  16T  model  diffuser  configuration.  The  full  potential  of  the  present 
formulation  should  be  further  validated  for  the  Mach  number  effect  when  additional  data 
for  the  16T  model  diffuser  are  obtained.  The  extension  to  three-dimensional  flows  also 
should  be  carried  out  with  a vorticity-velocity  formulation  because  most  flow  problems  are 
three-dimensional  in  nature. 
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NOMENCLATURE 

A Variable,  = V2k  y/v 

a,  b Constants 

C Coefficient  in  the  eddy  viscosity  model 

Cp  Pressure  coefficient 

H Step  size 

Aht  First  step  size 

j Integer 

K Coordinate  transformation  constant 

\ 

k Turbulent  kinetic  energy  (TKE) 

M Mach  number 

PI,  P2  Vector  path  function 
p Pressure 

Pexii,  w Exit  wall  pressure 

R Gas  constant 

r Radial  coordinate 

Sp  Pressure  source  term 

T Temperature 

u,  v Axial  and  radial  velocity  components 

Au  Wake  component  in  velocity  profile 

uc  Reference  velocity 

Uexit.c  Exit  reference  velocity 

ui(C  Inlet  reference  velocity 
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Unican  Mean  velocity 

v*  Friction  velocity 

x Axial  coordinate 

y Normal  coordinate 

Ay  Step  size  in  y 

Ay!  Boundary-layer  thickness 

6 2-D  flow,  = 0;  axisymmetric  flow,  = 1 

e TKE  dissipation 

H Molecular  viscosity 

fa  Eddy  viscosity 

ut  Eddy  viscosity  divided  by  density  Ot/p) 

p Density 

4/  Stream  function 

fl  Vorticity 
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